#------------------------------------------------------------------------------#
######################## Regressions Models ############################
#------------------------------------------------------------------------------#

rm(list = ls())

# Setting directory: adjust accordingly to 

setwd()

#------------------------------------------------------------------------------#
### Loading and installing necessary packages ##################################
#------------------------------------------------------------------------------#

# Installing required packages 

packages <- c('tidyverse', 'readxl', 'stringr', 
              'lmtest', 'sandwich', 'stargazer',
              'haven', 'survey', 'GDAtools',
              'weights','knitr','broom','gtools',
              'huxtable', 'compareGroups', 'lfe',
              'kableExtra', 'interplot', 'gridExtra',
              'miceadds', 'ggforce', 'httr', 'gdata',
              'rmarkdown', 'fixest', 'naniar', 'psych',
              'GPArotation', 'jtools', 'formattable',
              'openxlsx', 'fastDummies', 'mfx', 'margins',
              'corrplot', 'grid', 'gridExtra', 'broom.mixed',
              'erer', 'fixest', 'ggstance', 'lme4', 'lmerTest',
              'nlme', 'psych', 'spatstat', 'marginaleffects', 'texreg',
              'ggstance')

# Checking if is installed (and install if not)

packages.check <- lapply(
  packages,
  FUN = function(x) {
    if (!require(x, character.only = TRUE)) {
      install.packages(x, dependencies = TRUE)
      library(x, character.only = TRUE)
    }
  }
)


rm(packages, packages.check)

#------------------------------------------------------------------------------#
### Loading data base  #########################################################
#------------------------------------------------------------------------------#

Base <- readRDS(file = "data/Base.rds")

#------------------------------------------------------------------------------#
### Renaming Variables  ########################################################
#------------------------------------------------------------------------------#

Base <- Base %>%
  mutate(
    bf = `Belief climate change`,
    cau = `Causes of climate change`,
    imp = `Impacts of climate change`,
    pi_lf = `Political Ideology (Left-Right)`,
    pi_cs = `Political Ideology (Conservatism-Progressive)`, 
    sk = `Subjective knowledge`,
    ok = `Objective knowledge`,
    sc = `Scientific consensus`,
    ts = `Trust in scientists`, 
    nep = `The New Ecological Paradigm (NEP)`,
    ii = `Individualism index`,
    ei = `Egalitarianism index`,
    pe = `Personal Experience (extreme weather events)`)

#------------------------------------------------------------------------------#
### Creating Dictionary Variables  #############################################
#------------------------------------------------------------------------------#

mydict <- c("pi_lf" = "Political ideology: Left",  
            "pi_cs" = "Political ideology: Progressive",
            "sk" = "Subjective knowledge", 
            "ok" = "Objective knowledge",  
            "sc" = "Scientific consensus",
            "ts" = "Trust in scientists",
            "nep" = "The New Ecological Paradigm (NEP)",
            "ii" = "Individualism worldview",
            "ei" = "Egalitarianism worldview",
            "pe" = "Personal experience (extreme weather events)",
            "FemaleYes" = "Female", 
            "Escolaridade_3Highschoolorequivalent" = "Education: High school or equivalent", 
            "Escolaridade_3Undergraduateormore" = "Education: Undergraduate or more",
            "ReligionCatholic" = "Religion: Catholic", 
            "ReligionEvangelicalPentecostalorotherevangelical" = "Religion: Evangelical Pentecostal or other evangelical", 
            "ReligionEvangelicalTraditional" = "Religion: Evangelical Traditional",
            "ReligionOthers/NoRelig." = "Religion: Others/No Relig.",
            "Income1-2minimumwages" = "Income: 1 - 2 minimum wages", 
            "Income2-3minimumwages" = "Income: 2 - 3 minimum wages",
            "Income3-5minimumwages" = "Income: 3 - 5 minimum wages", 
            "Income5-10minimumwages" = "Income: 5 - 10 minimum wages",
            "Income10minimumwagesormore" = "Income: 10 minimum wages or more", 
            "IncomeDonotknow/Prefertonotanswer" = "Income: Do not know/ Prefer to not answer", 
            "age" = "Age (Years)", 
            "Color4" = "Race: Black",
            "(Intercept)" = "Const")

#------------------------------------------------------------------------------#
################ Regressions, Robustness and Tables ############################
#------------------------------------------------------------------------------#

# Setting baseline for socio-demographic variables
Base <- Base %>%
  mutate(Income = as.factor(Income),
         Income = factor(Income, levels = c("0 - 1 minimum wages", "1 - 2 minimum wages",
                                            "2 - 3 minimum wages", "3 - 5 minimum wages",
                                            "5 - 10 minimum wages", "10 minimum wages or more",
                                            "Do not know/ Prefer to not answer")),
         Income = relevel(Income, ref = "0 - 1 minimum wages"),
         Escolaridade_3 = as.factor(Escolaridade_3),
         Escolaridade_3 = relevel(Escolaridade_3, "Elementary (Primary) or less"),
         Religion = as.factor(Religion),
         Religion = relevel(Religion, "Atheist"),
         Female = as.factor(Female),
         Female = relevel(Female, "No"))

#------------------------------------------------------------------------------#
### (I) Existence of Climate Change (OLS) #########################################
#------------------------------------------------------------------------------#

regs <- list()  
for (i in seq(1,7)){ 
  regs[[i]] <- feols(bf ~ pi_lf + pi_cs + sk + ok + sc + ts + nep + ii + ei + pe + 
                       Female + Escolaridade_3 + Religion + Income + age + Color4, 
                     data=Base[Base$Pais==i,], weights = ~ ponde2, vcov = "hetero")
}

regs_total_I <- feols(bf ~ pi_lf + pi_cs + sk + ok + sc + ts + nep + ii + ei + pe + 
                        Female + Escolaridade_3 + Religion + Income + age + Color4, 
                      data = Base, weights = ~ ponde2, vcov = "hetero")

etable(regs[2],regs[1],regs[3],regs[4],regs[5],regs[6],regs[7],regs_total_I,
       title = "OLS Results - Existence of Climate Change",
       digits = 3,
       dict = mydict,
       signif.code = c(`***` = 0.01, `**` = 0.05, `*` = 0.1),
       file = "tables/ols_I.tex")

#------------------------------------------------------------------------------#
### (I) Existence of Climate Change (Ordinal Logit) ###############################
#------------------------------------------------------------------------------#

Base_ol <- Base %>% mutate(bf = factor(bf, ordered = TRUE, levels = c(0,1,2,3,4,5,6,7,8)))

arg <- polr(bf ~ pi_lf + pi_cs + sk + ok + sc + ts + nep + ii + ei + pe + 
              Female + Escolaridade_3 + Religion + Income + age + Color4, 
            data=Base_ol[Base_ol$Pais==2,],weights=ponde2, Hess = TRUE, start= c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
                                                                           0,0,0,0,0,1,2,3,4,5,6,7,8))
br <- polr(bf ~ pi_lf + pi_cs + sk + ok + sc + ts + nep + ii + ei + pe + 
             Female + Escolaridade_3 + Religion + Income + age + Color4, 
           data=Base_ol[Base_ol$Pais==1,],weights=ponde2, Hess = TRUE, start= c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
                                                                          0,0,0,0,0,1,2,3,4,5,6,7,8))
ch <- polr(bf ~ pi_lf + pi_cs + sk + ok + sc + ts + nep + ii + ei + pe + 
             Female + Escolaridade_3 + Religion + Income + age + Color4, 
           data=Base_ol[Base_ol$Pais==3,],weights=ponde2, Hess = TRUE, start= c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
                                                                          0,0,0,0,0,1,2,3,4,5,6,7,8))
co <- polr(bf ~ pi_lf + pi_cs + sk + ok + sc + ts + nep + ii + ei + pe + 
             Female + Escolaridade_3 + Religion + Income + age + Color4, 
           data=Base_ol[Base_ol$Pais==4,],weights=ponde2, Hess = TRUE, start= c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
                                                                          0,0,0,0,0,1,2,3,4,5,6,7,8))
eq <- polr(bf ~ pi_lf + pi_cs + sk + ok + sc + ts + nep + ii + ei + pe + 
             Female + Escolaridade_3 + Religion + Income + age + Color4, 
           data=Base_ol[Base_ol$Pais==5,],weights=ponde2, Hess = TRUE, start= c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
                                                                          0,0,0,0,0,1,2,3,4,5,6,7,8))
mx <- polr(bf ~ pi_lf + pi_cs + sk + ok + sc + ts + nep + ii + ei + pe + 
             Female + Escolaridade_3 + Religion + Income + age + Color4, 
           data=Base_ol[Base_ol$Pais==6,],weights=ponde2, Hess = TRUE, start= c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
                                                                          0,0,0,0,0,1,2,3,4,5,6,7,8))
pr <- polr(bf ~ pi_lf + pi_cs + sk + ok + sc + ts + nep + ii + ei + pe + 
             Female + Escolaridade_3 + Religion + Income + age + Color4, 
           data=Base_ol[Base_ol$Pais==7,],weights=ponde2, Hess = TRUE, start= c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
                                                                          0,0,0,0,0,1,2,3,4,5,6,7,8))

regs_total_logit_I <- polr(bf ~ pi_lf + pi_cs + sk + ok + sc + ts + nep + ii + ei + pe + 
                             Female + Escolaridade_3 + Religion + Income + age + Color4, 
                           data=Base_ol, weights=ponde2, Hess = TRUE, start= c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
                                                                            0,0,0,0,0,1,2,3,4,5,6,7,8))

texreg(list(arg, br, ch, co, eq, mx, pr, regs_total_logit_I),
       custom.coef.names = names)

#------------------------------------------------------------------------------#
### (II) Causes of Climate Change (OLS) ########################################
#------------------------------------------------------------------------------#

regs <- list()  
for (i in seq(1,7)){ 
  regs[[i]] <- feols(cau ~ pi_lf + pi_cs + sk + ok + sc + ts + nep + ii + ei + pe + 
                       Female + Escolaridade_3 + Religion + Income + age + Color4, 
                     data=Base[Base$Pais==i,], weights = ~ ponde2, vcov = "hetero")
}

regs_total_II <- feols(cau ~ pi_lf + pi_cs + sk + ok + sc + ts + nep + ii + ei + pe + 
                         Female + Escolaridade_3 + Religion + Income + age + Color4, 
                       data = Base, weights = ~ ponde2, vcov = "hetero")


etable(regs[2],regs[1],regs[3],regs[4],regs[5],regs[6],regs[7],regs_total_II,
       title = "OLS Results - Causes of Climate Change",
       digits = 3,
       dict = mydict,
       signifCode = c(`***` = 0.01, `**` = 0.05, `*` = 0.1),
       file = "tables/ols_II.tex")

#------------------------------------------------------------------------------#
### (II) Causes of Climate Change (Logit) ######################################
#------------------------------------------------------------------------------#

regs <- list()  
for (i in seq(1,7)){ 
  regs[[i]] <- feglm(cau ~ pi_lf + pi_cs + sk + ok + sc + ts + nep + ii + ei + pe + 
                       Female + Escolaridade_3 + Religion + Income + age + Color4, 
                     data=Base[Base$Pais==i,], family = "logit", weights = ~ ponde2, vcov = "hetero")
}

regs_total_logit_II <- feglm(cau ~ pi_lf + pi_cs + sk + ok + sc + ts + nep + ii + ei + pe + 
                               Female + Escolaridade_3 + Religion + Income + age + Color4, 
                             data = Base, family = "logit", weights = ~ ponde2, vcov = "hetero")


etable(regs[2],regs[1],regs[3],regs[4],regs[5],regs[6],regs[7],regs_total_logit_II,
       title = "Logit Results - Causes of Climate Change",
       digits = 3,
       dict = mydict,
       signifCode = c(`***` = 0.01, `**` = 0.05, `*` = 0.1),
       file = "tables/logit_II.tex")

#------------------------------------------------------------------------------#
### (III) Consequences of Climate Change (OLS) ######################################
#------------------------------------------------------------------------------#

regs <- list()  
for (i in seq(1,7)){ 
  regs[[i]] <- feols(imp ~ pi_lf + pi_cs + sk + ok + sc + ts + nep + ii + ei + pe + 
                       Female + Escolaridade_3 + Religion + Income + age + Color4, 
                     data=Base[Base$Pais==i,], weights = ~ ponde2, vcov = "hetero")
}

regs_total_III <- feols(imp ~ pi_lf + pi_cs + sk + ok + sc + ts + nep + ii + ei + pe + 
                          Female + Escolaridade_3 + Religion + Income + age + Color4, 
                        data = Base, weights = ~ ponde2, vcov = "hetero")


etable(regs[2],regs[1],regs[3],regs[4],regs[5],regs[6],regs[7],regs_total_III,
       title = "OLS Results - Impacts of Climate Change",
       digits = 3,
       dict = mydict,
       signifCode = c(`***` = 0.01, `**` = 0.05, `*` = 0.1),
       file = "tables/ols_III.tex")

#------------------------------------------------------------------------------#
### (III) Consequences of Climate Change (Logit) ####################################
#------------------------------------------------------------------------------#

regs <- list()  
for (i in seq(1,7)){ 
  regs[[i]] <- feglm(imp ~ pi_lf + pi_cs + sk + ok + sc + ts + nep + ii + ei + pe + 
                       Female + Escolaridade_3 + Religion + Income + age + Color4, 
                     data=Base[Base$Pais==i,], family = "logit", weights = ~ ponde2, vcov = "hetero")
}

regs_total_logit_III <- feglm(imp ~ pi_lf + pi_cs + sk + ok + sc + ts + nep + ii + ei + pe + 
                                Female + Escolaridade_3 + Religion + Income + age + Color4, 
                              data = Base, family = "logit", weights = ~ ponde2, vcov = "hetero")


etable(regs[2],regs[1],regs[3],regs[4],regs[5],regs[6],regs[7],regs_total_logit_III,
       title = "Logit Results - Impacts of Climate Change",
       digits = 3,
       dict = mydict,
       signifCode = c(`***` = 0.01, `**` = 0.05, `*` = 0.1),
       file = "tables/logit_III.tex")

#------------------------------------------------------------------------------#
### (IV) Belief in Climate Change Index (OLS) ######################################
#------------------------------------------------------------------------------#

regs <- list()  
for (i in seq(1,7)){ 
  regs[[i]] <- feols(index_std ~ pi_lf + pi_cs + sk + ok + sc + ts + nep + ii + ei + pe + 
                       Female + Escolaridade_3 + Religion + Income + age + Color4, 
                     data=Base[Base$Pais==i,], weights = ~ ponde2, vcov = "hetero")
}

regs_total_IV <- feols(index_std ~ pi_lf + pi_cs + sk + ok + sc + ts + nep + ii + ei + pe + 
                          Female + Escolaridade_3 + Religion + Income + age + Color4, 
                        data = Base, weights = ~ ponde2, vcov = "hetero")


etable(regs[2],regs[1],regs[3],regs[4],regs[5],regs[6],regs[7],regs_total_IV,
       title = "OLS Results - Belief in Climate Change (Index)",
       digits = 3,
       dict = mydict,
       signifCode = c(`***` = 0.01, `**` = 0.05, `*` = 0.1),
       file = "tables/ols_IV.tex")

#------------------------------------------------------------------------------#
### (I-IV) Multilevel Model ###################################################
#------------------------------------------------------------------------------#

# (I) Belief in Climate Change 

regs_total_I_ml <- lmer(bf ~  pi_lf + pi_cs + sk + ok + sc + ts + nep + ii + ei + pe + 
                          Female + Escolaridade_3 + Religion + Income + age + Color4 + (1|Pais), data=Base, weights=ponde2)
summary(regs_total_I_ml)

# (II) Causes of Climate Change 

regs_total_II_ml <- lmer(cau ~  pi_lf + pi_cs + sk + ok + sc + ts + nep + ii + ei + pe + 
                           Female + Escolaridade_3 + Religion + Income + age + Color4 + (1|Pais), data=Base,weights=ponde2)
summary(regs_total_II_ml)

# (III) Impacts of Climate Change 

regs_total_III_ml <- lmer(imp ~  pi_lf + pi_cs + sk + ok + sc + ts + nep + ii + ei + pe + 
                            Female + Escolaridade_3 + Religion + Income + age + Color4 + (1|Pais), data=Base,weights=ponde2)
summary(regs_total_III_ml)

# Table Models: I-III

texreg(list(regs_total_I_ml, regs_total_II_ml, regs_total_III_ml),
       custom.coef.names = names_c)

#------------------------------------------------------------------------------#
### (I-IV) Fixed Effects Models ###############################################
#------------------------------------------------------------------------------#

# (I) Existence of Climate Change 
regs_total_I_fe <- feols(bf ~ pi_lf + pi_cs + sk + ok + sc + ts + nep + ii + ei + pe + 
                           Female + Escolaridade_3 + Religion + Income + age + Color4 | Pais, 
                         data = Base, weights = ~ ponde2, vcov = "hetero")

# (II) Causes of Climate Change 
regs_total_II_fe <- feols(cau ~ pi_lf + pi_cs + sk + ok + sc + ts + nep + ii + ei + pe + 
                            Female + Escolaridade_3 + Religion + Income + age + Color4 | Pais, 
                          data = Base, weights = ~ ponde2, vcov = "hetero")

# (III) Impacts of Climate Change 
regs_total_III_fe <- feols(imp ~ pi_lf + pi_cs + sk + ok + sc + ts + nep + ii + ei + pe + 
                             Female + Escolaridade_3 + Religion + Income + age + Color4 | Pais, 
                           data = Base, weights = ~ ponde2, vcov = "hetero")

# (IV) Belief in Climate Change 
regs_total_IV_fe <- feols(index_std ~ pi_lf + pi_cs + sk + ok + sc + ts + nep + ii + ei + pe + 
                             Female + Escolaridade_3 + Religion + Income + age + Color4 | Pais, 
                           data = Base, weights = ~ ponde2, vcov = "hetero")

etable(regs_total_I_fe, regs_total_II_fe, regs_total_III_fe, regs_total_IV_fe)

etable(regs_total_I_fe, regs_total_II_fe, regs_total_III_fe, regs_total_IV_fe,
       title = "FE Results",
       digits = 3,
       dict = mydict,
       signif.code = c(`***` = 0.01, `**` = 0.05, `*` = 0.1),
       file = "tables/fe_I_II_III_IV.tex")

#------------------------------------------------------------------------------#
################ Plotting Results ##############################################
#------------------------------------------------------------------------------#

# Naming variables 

mydict_plot <- c("Political ideology: Left" = "pi_lf",
            "Political ideology: Progressive" ="pi_cs",
            "Subjective knowledge" = "sk",
            "Objective knowledge" = "ok",
            "Scientific consensus" = "sc",
            "Trust in scientists" = "ts",
            "The New Ecological Paradigm (NEP)" = "nep",
            "Individualism worldview" = "ii",
            "Egalitarianism worldview" =  "ei",
            "Personal experience (extreme weather events)" = "pe",
            "Female" = "FemaleYes", 
            "Education: High school or equivalent" = "Escolaridade_3Highschoolorequivalent", 
            "Education: Undergraduate or more" = "Escolaridade_3Undergraduateormore",
            "Religion: Catholic" = "ReligionCatholic", 
            "Religion: Evangelical Pentecostal or other evangelical" = "ReligionEvangelicalPentecostalorotherevangelical", 
            "Religion: Evangelical Traditional" = "ReligionEvangelicalTraditional",
            "Religion: Others/No Relig." = "ReligionOthers/NoRelig.",
            "Income: 1 - 2 minimum wages" = "Income1-2minimumwages", 
            "Income: 2 - 3 minimum wages" = "Income2-3minimumwages",
            "Income: 3 - 5 minimum wages" = "Income3-5minimumwages", 
            "Income: 5 - 10 minimum wages" = "Income5-10minimumwages",
            "Income: 10 minimum wages or more" = "Income10minimumwagesormore", 
            "Income: Do not know/ Prefer to not answer" = "IncomeDonotknow/Prefertonotanswer", 
            "Age (Years)" = "age", 
            "Race: Black" = "Color4")

names <- c("Political ideology: Left",  "Political ideology: Progressive",
             "Subjective knowledge", "Objective knowledge",  "Scientific consensus","Trust in scientists",
             "The New Ecological Paradigm (NEP)", "Individualism worldview","Egalitarianism worldview",
             "Personal experience (extreme weather events)",
             "Female", "Education: High school or equivalent", "Education: Undergraduate or more",
             "Religion: Catholic", "Religion: Evangelical Pentecostal or other evangelical", "Religion: Evangelical Traditional",
             "Religion: Others/No Relig.",
             "Income: 1 - 2 minimum wages", "Income: 2 - 3 minimum wages",
             "Income: 3 - 5 minimum wages", "Income: 5 - 10 minimum wages",
             "Income: 10 minimum wages or more", "Income: Do not know/ Prefer to not answer", "Age (Years)", "Race: Black", "0", "1", "2", "3", "4", "5", "6", "7")

names_c <- c("Constant", "Political ideology: Left",  "Political ideology: Progressive",
             "Subjective knowledge", "Objective knowledge",  "Scientific consensus","Trust in scientists",
             "The New Ecological Paradigm (NEP)", "Individualism worldview","Egalitarianism worldview",
             "Personal experience (extreme weather events)",
             "Female", "Education: High school or equivalent", "Education: Undergraduate or more",
             "Religion: Catholic", "Religion: Evangelical Pentecostal or other evangelical", "Religion: Evangelical Traditional",
             "Religion: Others/No Relig.",
             "Income: 1 - 2 minimum wages", "Income: 2 - 3 minimum wages",
             "Income: 3 - 5 minimum wages", "Income: 5 - 10 minimum wages",
             "Income: 10 minimum wages or more", "Income: Do not know/ Prefer to not answer", "Age (Years)", "Race: Black")

indep <- c("Subjective knowledge" = "sk",
            "Objective knowledge" = "ok",
            "Scientific consensus" = "sc",
            "Trust in scientists" = "ts",
            "The New Ecological Paradigm (NEP)" = "nep",
            "Individualism worldview" = "ii",
            "Egalitarianism worldview" =  "ei",
            "Personal experience (extreme weather events)" = "pe")

indep_sub <- c("Subjective knowledge",
           "Objective knowledge",
           "Scientific consensus",
           "Trust in scientists",
           "The New Ecological Paradigm (NEP)",
           "Individualism worldview",
           "Egalitarianism worldview",
           "Personal experience (extreme weather events)")

socio_dem <- c("Political ideology: Left" = "pi_lf",
               "Political ideology: Progressive" ="pi_cs",
               "Female" = "FemaleYes", 
               "Education: High school or equivalent" = "Escolaridade_3High school or equivalent", 
               "Education: Undergraduate or more" = "Escolaridade_3Undergraduate or more",
               "Religion: Catholic" = "ReligionCatholic", 
               "Religion: Evangelical Pentecostal or other evangelical" = "ReligionEvangelical Pentecostal or other evangelical", 
               "Religion: Evangelical Traditional" = "ReligionEvangelical Traditional",
               "Religion: Others/No Relig." = "ReligionOthers/No Relig.",
               "Income: 1 - 2 minimum wages" = "Income1 - 2 minimum wages", 
               "Income: 2 - 3 minimum wages" = "Income2 - 3 minimum wages",
               "Income: 3 - 5 minimum wages" = "Income3 - 5 minimum wages", 
               "Income: 5 - 10 minimum wages" = "Income5 - 10 minimum wages",
               "Income: 10 minimum wages or more" = "Income10 minimum wages or more", 
               "Income: Do not know/ Prefer to not answer" = "IncomeDo not know/ Prefer to not answer", 
               "Age (Years)" = "age", 
               "Race: Black" = "Color4")

socio_dem_sub <- c("Political ideology: Left",
                   "Political ideology: Progressive",
                   "Female", 
                   "Education: High school or equivalent", 
                   "Education: Undergraduate or more",
                   "Religion: Catholic", 
                   "Religion: Evangelical Pentecostal or other evangelical", 
                   "Religion: Evangelical Traditional",
                   "Religion: Others/No Relig.",
                   "Income: 1 - 2 minimum wages", 
                   "Income: 2 - 3 minimum wages",
                   "Income: 3 - 5 minimum wages", 
                   "Income: 5 - 10 minimum wages",
                   "Income: 10 minimum wages or more", 
                   "Income: Do not know/ Prefer to not answer", 
                   "Age (Years)", 
                   "Race: Black")

#------------------------------------------------------------------------------#
### Plot: (I) Existence of Climate Change (OLS) ###################################
#------------------------------------------------------------------------------#

label_regs_total_I <- tidy(regs_total_I, conf.int = TRUE, conf.level = 0.95) %>%
  mutate_if(is.numeric, round, digits = 2) %>%
  mutate(coef_ci = paste0(format(round(estimate, 2), nsmall = 2), " ", "[" , 
                          format(round(conf.low, 2), nsmall = 2), ",",
                          format(round(conf.high, 2), nsmall = 2), "]")) %>%
  mutate(term = names_c)

# Independent Variables
indep_I <- plot_coefs(regs_total_I, scale = TRUE,
                      coefs = indep, robust = TRUE) +
  xlim(c(-0.45, 0.68)) +
  geom_text(inherit.aes = FALSE, 
            data = subset(label_regs_total_I, term %in% indep_sub),
            aes(x = estimate, y = term, label = coef_ci), 
            size = 3.2, 
            vjust = -1) +
  ggtitle("\n Psychological variables") +
  theme_bw() +
  theme(axis.text.y = element_text(color="black", size = 11),
        axis.text.x = element_text(color="black"),
        axis.title.y = element_blank(),
        axis.title.x = element_text(size = 9),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())

# Socio-economic Variables
socio_dem_I <- plot_coefs(regs_total_I, scale = TRUE,
                          coefs = socio_dem, robust = TRUE) +
  xlim(c(-0.45, 0.68)) +
  geom_text(inherit.aes = FALSE, 
            data = subset(label_regs_total_I, term %in% socio_dem_sub),
            aes(x = estimate, y = term, label = coef_ci), 
            size = 3.2, 
            vjust = -1) +
  ggtitle("Political ideology and \n Socio-Demographic variables") +
  theme_bw() +
  theme(axis.text.y = element_text(color="black", size = 11),
        axis.text.x = element_text(color="black"),
        axis.title.y = element_blank(),
        axis.title.x = element_text(size = 9),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())

# Joining both plots
grid.arrange(indep_I, socio_dem_I,
             nrow = 1,
             widths = c(1,1.08)) #Salvar como PDF A4 landscape width 12.5 e height 7 e add como pdf no overleaf

#------------------------------------------------------------------------------#
### Plot: (II) Causes of Climate Change (OLS) ##################################
#------------------------------------------------------------------------------#

label_regs_total_II <- tidy(regs_total_II, conf.int = TRUE, conf.level = 0.95) %>%
  mutate_if(is.numeric, round, digits = 2) %>%
  mutate(coef_ci = paste0(format(round(estimate, 2), nsmall = 2), " ", "[" , 
                          format(round(conf.low, 2), nsmall = 2), ",",
                          format(round(conf.high, 2), nsmall = 2), "]")) %>%
  mutate(term = names_c)

# Independent Variables
indep_II <- plot_coefs(regs_total_II, scale = TRUE,
                       coefs = indep, robust = TRUE) +
  xlim(c(-0.2, 0.15)) +
  geom_text(inherit.aes = FALSE, 
            data = subset(label_regs_total_II, term %in% indep_sub),
            aes(x = estimate, y = term, label = coef_ci), 
            size = 3.2, 
            vjust = -1) +
  ggtitle("\n Psychological variables") +
  theme_bw() +
  theme(axis.text.y = element_text(color="black", size = 11),
        axis.text.x = element_text(color="black"),
        axis.title.y = element_blank(),
        axis.title.x = element_text(size = 9),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())

# Socio-economic Variables
socio_dem_II <- plot_coefs(regs_total_II, scale = TRUE,
                           coefs = socio_dem, robust = TRUE) +
  xlim(c(-0.2, 0.15)) +
  geom_text(inherit.aes = FALSE, 
            data = subset(label_regs_total_II, term %in% socio_dem_sub),
            aes(x = estimate, y = term, label = coef_ci), 
            size = 3.2, 
            vjust = -1) +
  ggtitle("Political ideology and \n Socio-Demographic variables") +
  theme_bw() +
  theme(axis.text.y = element_text(color="black", size = 11),
        axis.text.x = element_text(color="black"),
        axis.title.y = element_blank(),
        axis.title.x = element_text(size = 9),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())

# Joining both plots
grid.arrange(indep_II, socio_dem_II,
             nrow = 1,
             widths = c(1,1.05)) #Salvar como PDF A4 landscape width 12.5 e height 7 e add como pdf no overleaf

#------------------------------------------------------------------------------#
### Plot: (III) Consequences of Climate Change (OLS) ################################
#------------------------------------------------------------------------------#

label_regs_total_III <- tidy(regs_total_III, conf.int = TRUE, conf.level = 0.95) %>%
  mutate_if(is.numeric, round, digits = 2) %>%
  mutate(coef_ci = paste0(format(round(estimate, 2), nsmall = 2), " ", "[" , 
                          format(round(conf.low, 2), nsmall = 2), ",",
                          format(round(conf.high, 2), nsmall = 2), "]")) %>%
  mutate(term = names_c)

# Independent Variables
indep_III <- plot_coefs(regs_total_III, scale = TRUE,
                        coefs = indep, robust = TRUE) +
  xlim(c(-0.3, 0.25)) +
  geom_text(inherit.aes = FALSE, 
            data = subset(label_regs_total_III, term %in% indep_sub),
            aes(x = estimate, y = term, label = coef_ci), 
            size = 3.2, 
            vjust = -1) +
  ggtitle("\n Psychological variables") +
  theme_bw() +
  theme(axis.text.y = element_text(color="black", size = 11),
        axis.text.x = element_text(color="black"),
        axis.title.y = element_blank(),
        axis.title.x = element_text(size = 9),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())

# Socio-economic Variables
socio_dem_III <- plot_coefs(regs_total_III, scale = TRUE,
                            coefs = socio_dem, robust = TRUE) +
  xlim(c(-0.3, 0.25)) +
  geom_text(inherit.aes = FALSE, 
            data = subset(label_regs_total_III, term %in% socio_dem_sub),
            aes(x = estimate, y = term, label = coef_ci), 
            size = 3.2, 
            vjust = -1) +
  ggtitle("Political ideology and \n Socio-Demographic variables") +
  theme_bw() +
  theme(axis.text.y = element_text(color="black", size = 11),
        axis.text.x = element_text(color="black"),
        axis.title.y = element_blank(),
        axis.title.x = element_text(size = 9),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())

# Joining both plots
grid.arrange(indep_III, socio_dem_III,
             nrow = 1, 
             widths = c(1,1.05)) #Salvar como PDF A4 landscape width 12.5 e height 7 e add como pdf no overleaf

#------------------------------------------------------------------------------#
### Plot: (II-III) (OLS) #######################################################
#------------------------------------------------------------------------------#

# Independent Variables without legend 
indep_II_III <- plot_coefs(regs_total_II, regs_total_III,
                           model.names= c("Anthropogenic Climate Change",
                                          "Consequences of Climate Change"),
                           legend.title = "Outcomes",
                           coefs = indep,
                           robust = TRUE) +
  xlim(c(-0.3, 0.25)) +
  ggtitle("\n Psychological variables") +
  theme_bw() +
  theme(axis.text.y = element_text(color="black", size = 11),
        axis.text.x = element_text(color="black"),
        axis.title.y = element_blank(),
        axis.title.x = element_text(size = 9),
        legend.position = "none",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())

# Socio-economic Variables without legend 
socio_dem_II_III <- plot_coefs(regs_total_II, regs_total_III,
                               model.names= c("Anthropogenic Climate Change",
                                              "Consequences of Climate Change"),
                               legend.title = "Outcomes",
                               coefs = socio_dem,
                               robust = TRUE) +
  xlim(c(-0.3, 0.25)) +
  ggtitle("Political ideology and \n Socio-Demographic variables") +
  theme_bw() +
  theme(axis.text.y = element_text(color="black", size = 11),
        axis.text.x = element_text(color="black"),
        axis.title.y = element_blank(),
        axis.title.x = element_text(size = 9),
        legend.position = "none",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())

# Joining both plots 
combined_plot <- grid.arrange(indep_II_III, socio_dem_II_III,
                              nrow = 1,
                              widths = c(1,1.05))

# Independent Variables with legend
indep_II_III_legend <- plot_coefs(regs_total_II, regs_total_III,
                                  model.names= c("Anthropogenic Climate Change",
                                                 "Consequences of Climate Change"),
                                  legend.title = "Outcome",
                                  coefs = indep,
                                  robust = TRUE) +
  xlim(c(-0.25, 0.25)) +
  ggtitle("\n Psychological variables") +
  theme_bw() +
  theme(axis.text.y = element_text(color="black", size = 11),
        axis.text.x = element_text(color="black"),
        axis.title.y = element_blank(),
        axis.title.x = element_text(size = 9),
        legend.position = "bottom")

# Function to extract legend
get_only_legend <- function(plot) {
  plot_table <- ggplot_gtable(ggplot_build(plot))
  legend_plot <- which(sapply(plot_table$grobs, function(x) x$name) == "guide-box")
  legend <- plot_table$grobs[[legend_plot]]
  return(legend)
}

# Extract legend from Independent Variables with legend using above function
legend <- get_only_legend(indep_II_III_legend)

# Joining both plots wit legends
grid.arrange(combined_plot,
             legend,
             nrow = 2,
             heights = c(10, 1)) #Salvar como PDF A4 landscape width 12.5 e height 7 e add como pdf no overleaf


#------------------------------------------------------------------------------#
### Plot: (IV) Belief in Climate Change Index (OLS) ################################
#------------------------------------------------------------------------------#

label_regs_total_IV <- tidy(regs_total_IV, conf.int = TRUE, conf.level = 0.95) %>%
  mutate_if(is.numeric, round, digits = 2) %>%
  mutate(coef_ci = paste0(format(round(estimate, 2), nsmall = 2), " ", "[" , 
                          format(round(conf.low, 2), nsmall = 2), ",",
                          format(round(conf.high, 2), nsmall = 2), "]")) %>%
  mutate(term = names_c)

# Independent Variables
indep_IV <- plot_coefs(regs_total_IV, scale = TRUE,
                      coefs = indep, robust = TRUE) +
  xlim(c(-0.45, 0.68)) +
  geom_text(inherit.aes = FALSE, 
            data = subset(label_regs_total_IV, term %in% indep_sub),
            aes(x = estimate, y = term, label = coef_ci), 
            size = 3.2, 
            vjust = -1) +
  ggtitle("\n Psychological variables") +
  theme_bw() +
  theme(axis.text.y = element_text(color="black", size = 11),
        axis.text.x = element_text(color="black"),
        axis.title.y = element_blank(),
        axis.title.x = element_text(size = 9),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())

# Socio-economic Variables
socio_dem_IV <- plot_coefs(regs_total_IV, scale = TRUE,
                          coefs = socio_dem, robust = TRUE) +
  xlim(c(-0.65, 0.68)) +
  geom_text(inherit.aes = FALSE, 
            data = subset(label_regs_total_IV, term %in% socio_dem_sub),
            aes(x = estimate, y = term, label = coef_ci), 
            size = 3.2, 
            vjust = -1) +
  ggtitle("Political ideology and \n Socio-Demographic variables") +
  theme_bw() +
  theme(axis.text.y = element_text(color="black", size = 11),
        axis.text.x = element_text(color="black"),
        axis.title.y = element_blank(),
        axis.title.x = element_text(size = 9),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())

# Joining both plots
grid.arrange(indep_IV, socio_dem_IV,
             nrow = 1,
             widths = c(1,1.08)) #Salvar como PDF A4 landscape width 12.5 e height 7 e add como pdf no overleaf

#------------------------------------------------------------------------------#
### Plotting OLS per Country ###################################################
#------------------------------------------------------------------------------#

demog <- c("pi_lf", "pi_cs", "FemaleYes", "Escolaridade_3High school or equivalent", "Escolaridade_3Undergraduate or more", 
           "ReligionCatholic", "ReligionEvangelical Pentecostal or other evangelical", "ReligionEvangelical Traditional", 
           "ReligionOthers/No Relig.", "Income1 - 2 minimum wages", "Income2 - 3 minimum wages", 
           "Income3 - 5 minimum wages", "Income5 - 10 minimum wages", "Income10 minimum wages or more", 
           "IncomeDo not know/ Prefer to not answer", "age", "Color4")

indepen <- c("sk", "ok", "sc", "ts", "nep", "ii", "ei", "pe")

countries <- c("Brazil", "Argentina", "Chile", "Colombia", "Ecuador", "Mexico", "Peru")

#------------------------------------------------------------------------------#
### Plot (I) : Existence of Climate Change per Country (OLS) ######################
#------------------------------------------------------------------------------#

regI <- feols(bf ~ pi_lf + pi_cs + sk + ok + sc + ts + nep + ii + ei + pe +
              Female + Escolaridade_3 + Religion + Income + age + Color4, 
              data=Base, weights= ~ ponde2, split = ~ Pais, vcov = "hetero")

names(regI) <- countries

# Saving data frame with estimates and ci of model for each country

# Creating function to add CI to tidy function
tidy_ci <- function(x){
  tidy(x, conf.int = TRUE, conf.level = 0.95)
} 

regI_df <- map_dfr(regI, tidy_ci, .id = "Pais") %>%
  filter(term != "(Intercept)")

# Plotting results for Independent variables
indepen_I <- regI_df %>%
  mutate(term = factor(term, level=c("pe", "ei", "ii", "nep", "ts", "sc", "ok", "sk"))) %>%
  filter(term %in% indepen) %>%
  mutate(term = dplyr::recode(term,
                       "sk" = "Subjective knowledge", 
                       "ok" = "Objective knowledge", 
                       "sc" = "Scientific consensus",
                       "ts" = "Trust in scientists", 
                       "nep" = "The New Ecological Paradigm (NEP)", 
                       "ii" = "Individualism worldview", 
                       "ei" = "Egalitarianism worldview", 
                       "pe" = "Personal experience (extreme weather events)")) %>%
  ggplot(aes(x = estimate, y = term)) +
  geom_vline(xintercept = 0, colour = gray(1/2), lty = 2, size = 0.25) +
  geom_linerange(aes(x = estimate, 
                     xmin = conf.low,
                     xmax = conf.high),
                 lwd = 0.65) +
  geom_point(aes(x = estimate, 
                 y = term),
             shape = 21,
             size = 2,
             fill = "white") +
  xlab("Estimate") +
  ggtitle("Psychological variables") +
  theme_bw() +
  theme(axis.text.y = element_text(color="black", size = 10),
        axis.text.x = element_text(color="black"),
        axis.title.y = element_blank(),
        axis.title.x = element_text(size = 9),
        legend.position = "none",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()) +
  facet_grid(. ~ Pais)


# Plotting results for Socio-demographic variables
sociodem_I <- regI_df %>%
  mutate(term = factor(term, level=c("Color4", "age", "IncomeDo not know/ Prefer to not answer", "Income10 minimum wages or more", 
                                     "Income5 - 10 minimum wages","Income3 - 5 minimum wages", "Income2 - 3 minimum wages", "Income1 - 2 minimum wages", 
                                     "ReligionOthers/No Relig.", "ReligionEvangelical Traditional", "ReligionEvangelical Pentecostal or other evangelical", 
                                     "ReligionCatholic", "Escolaridade_3Undergraduate or more", "Escolaridade_3High school or equivalent", "FemaleYes", "pi_cs", "pi_lf"))) %>%
  filter(term %in% demog) %>%
  mutate(term = dplyr::recode(term,
                       "pi_lf" = "Political ideology: Left", 
                       "pi_cs" = "Political ideology: Progressive", 
                       "Color4" = "Race: Black", 
                       "age" = "Age (Years)", 
                       "IncomeDo not know/ Prefer to not answer" = "Income: Do not know / Prefer to not answer", 
                       "Income10 minimum wages or more" = "Income: 10 minimum wages or more", 
                       "Income5 - 10 minimum wages" = "Income: 5 - 10 minimum wages",
                       "Income3 - 5 minimum wages" = "Income: 3 - 5 minimum wages",
                       "Income2 - 3 minimum wages" = "Income: 2 - 3 minimum wages", 
                       "Income1 - 2 minimum wages" = "Income: 1 - 2 minimum wages", 
                       "ReligionOthers/No Relig." = "Religion: Others/No Relig.", 
                       "ReligionEvangelical Traditional" = "Religion: Evangelical Traditional", 
                       "ReligionEvangelical Pentecostal or other evangelical" = "Religion: Evangelical Pentecostal or other", 
                       "ReligionCatholic" = "Religion: Catholic",
                       "Escolaridade_3Undergraduate or more" = "Education: Undergraduate or more", 
                       "Escolaridade_3High school or equivalent" = "Education: High school or equivalent", 
                       "FemaleYes" = "Female")) %>%
  ggplot(aes(x = estimate, y = term)) +
  geom_vline(xintercept = 0, colour = gray(1/2), lty = 2, size = 0.25) +
  geom_linerange(aes(x = estimate, 
                     xmin = conf.low,
                     xmax = conf.high),
                 lwd = 0.65) +
  geom_point(aes(x = estimate, 
                 y = term),
             shape = 21,
             size = 2,
             fill = "white") +
  xlab("Estimate") +
  ggtitle("Political ideology and \n Socio-Demographic variables") +
  theme_bw() +
  theme(axis.text.y = element_text(color="black", size = 10),
        axis.text.x = element_text(color="black"),
        axis.title.y = element_blank(),
        axis.title.x = element_text(size = 9),
        legend.position = "none",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()) +
  facet_grid(. ~ Pais)

# Joining both plots
grid.arrange(indepen_I, sociodem_I,
             ncol = 1) # 8 x 12.5

#------------------------------------------------------------------------------#
### Plot (II) : Causes of Climate Change per Country (OLS) #####################
#------------------------------------------------------------------------------#

regII <- feols(cau ~ pi_lf + pi_cs + sk + ok + sc + ts + nep + ii + ei + pe +
                 Female + Escolaridade_3 + Religion + Income + age + Color4, 
               data=Base, weights= ~ ponde2, split = ~ Pais, vcov = "hetero")

names(regII) <- countries

# Saving data frame with estimates and ci of model for each country

# Creating function to add CI to tidy function
tidy_ci <- function(x){
  tidy(x, conf.int = TRUE, conf.level = 0.95)
} 

regII_df <- map_dfr(regII, tidy_ci, .id = "Pais") %>%
  filter(term != "(Intercept)")

# Plotting results for Independent variables
indepen_II <- regII_df %>%
  mutate(term = factor(term, level=c("pe", "ei", "ii", "nep", "ts", "sc", "ok", "sk"))) %>%
  filter(term %in% indepen) %>%
  mutate(term = dplyr::recode(term,
                       "sk" = "Subjective knowledge", 
                       "ok" = "Objective knowledge", 
                       "sc" = "Scientific consensus",
                       "ts" = "Trust in scientists", 
                       "nep" = "The New Ecological Paradigm (NEP)", 
                       "ii" = "Individualism worldview", 
                       "ei" = "Egalitarianism worldview", 
                       "pe" = "Personal experience (extreme weather events)")) %>%
  ggplot(aes(x = estimate, y = term)) +
  geom_vline(xintercept = 0, colour = gray(1/2), lty = 2, size = 0.25) +
  geom_linerange(aes(x = estimate, 
                     xmin = conf.low,
                     xmax = conf.high),
                 lwd = 0.65) +
  geom_point(aes(x = estimate, 
                 y = term),
             shape = 21,
             size = 2,
             fill = "white") +
  xlab("Estimate") +
  ggtitle("Psychological variables") +
  theme_bw() +
  theme(axis.text.y = element_text(color="black", size = 10),
        axis.text.x = element_text(color="black"),
        axis.title.y = element_blank(),
        axis.title.x = element_text(size = 9),
        legend.position = "none",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()) +
  facet_grid(. ~ Pais)


# Plotting results for Socio-demographic variables
sociodem_II <- regII_df %>%
  mutate(term = factor(term, level=c("Color4", "age", "IncomeDo not know/ Prefer to not answer", "Income10 minimum wages or more", 
                                     "Income5 - 10 minimum wages","Income3 - 5 minimum wages", "Income2 - 3 minimum wages", "Income1 - 2 minimum wages", 
                                     "ReligionOthers/No Relig.", "ReligionEvangelical Traditional", "ReligionEvangelical Pentecostal or other evangelical", 
                                     "ReligionCatholic", "Escolaridade_3Undergraduate or more", "Escolaridade_3High school or equivalent", "FemaleYes", "pi_cs", "pi_lf"))) %>%
  filter(term %in% demog) %>%
  mutate(term = dplyr::recode(term,
                       "pi_lf" = "Political ideology: Left", 
                       "pi_cs" = "Political ideology: Progressive", 
                       "Color4" = "Race: Black", 
                       "age" = "Age (Years)", 
                       "IncomeDo not know/ Prefer to not answer" = "    Income: Do not know / Prefer to not answer", 
                       "Income10 minimum wages or more" = "Income: 10 minimum wages or more", 
                       "Income5 - 10 minimum wages" = "Income: 5 - 10 minimum wages",
                       "Income3 - 5 minimum wages" = "Income: 3 - 5 minimum wages",
                       "Income2 - 3 minimum wages" = "Income: 2 - 3 minimum wages", 
                       "Income1 - 2 minimum wages" = "Income: 1 - 2 minimum wages", 
                       "ReligionOthers/No Relig." = "Religion: Others/No Relig.", 
                       "ReligionEvangelical Traditional" = "Religion: Evangelical Traditional", 
                       "ReligionEvangelical Pentecostal or other evangelical" = "Religion: Evangelical Pentecostal or other", 
                       "ReligionCatholic" = "Religion: Catholic",
                       "Escolaridade_3Undergraduate or more" = "Education: Undergraduate or more", 
                       "Escolaridade_3High school or equivalent" = "Education: High school or equivalent", 
                       "FemaleYes" = "Female")) %>%
  ggplot(aes(x = estimate, y = term)) +
  geom_vline(xintercept = 0, colour = gray(1/2), lty = 2, size = 0.25) +
  geom_linerange(aes(x = estimate, 
                     xmin = conf.low,
                     xmax = conf.high),
                 lwd = 0.65) +
  geom_point(aes(x = estimate, 
                 y = term),
             shape = 21,
             size = 2,
             fill = "white") +
  xlab("Estimate") +
  ggtitle("Political ideology and \n Socio-Demographic variables") +
  theme_bw() +
  theme(axis.text.y = element_text(color="black", size = 10),
        axis.text.x = element_text(color="black"),
        axis.title.y = element_blank(),
        axis.title.x = element_text(size = 9),
        legend.position = "none",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()) +
  facet_grid(. ~ Pais)

# Joining both plots
grid.arrange(indepen_II, sociodem_II,
             ncol = 1)

#------------------------------------------------------------------------------#
### Plot (III): Consequences of Climate Change per Country (OLS) ####################
#------------------------------------------------------------------------------#

regIII <- feols(imp ~ pi_lf + pi_cs + sk + ok + sc + ts + nep + ii + ei + pe +
                  Female + Escolaridade_3 + Religion + Income + age + Color4, 
                data=Base, weights= ~ ponde2, split = ~ Pais, vcov = "hetero")

names(regIII) <- countries

# Saving data frame with estimates and ci of model for each country

# Creating function to add CI to tidy function
tidy_ci <- function(x){
  tidy(x, conf.int = TRUE, conf.level = 0.95)
} 

regIII_df <- map_dfr(regIII, tidy_ci, .id = "Pais") %>%
  filter(term != "(Intercept)")

# Plotting results for Independent variables
indepen_III <- regIII_df %>%
  mutate(term = factor(term, level=c("pe", "ei", "ii", "nep", "ts", "sc", "ok", "sk"))) %>%
  filter(term %in% indepen) %>%
  mutate(term = dplyr::recode(term,
                       "sk" = "Subjective knowledge", 
                       "ok" = "Objective knowledge", 
                       "sc" = "Scientific consensus",
                       "ts" = "Trust in scientists", 
                       "nep" = "The New Ecological Paradigm (NEP)", 
                       "ii" = "Individualism worldview", 
                       "ei" = "Egalitarianism worldview", 
                       "pe" = "Personal experience (extreme weather events)")) %>%
  ggplot(aes(x = estimate, y = term)) +
  geom_vline(xintercept = 0, colour = gray(1/2), lty = 2, size = 0.25) +
  geom_linerange(aes(x = estimate, 
                     xmin = conf.low,
                     xmax = conf.high),
                 lwd = 0.65) +
  geom_point(aes(x = estimate, 
                 y = term),
             shape = 21,
             size = 2,
             fill = "white") +
  xlab("Estimate") +
  ggtitle("Psychological variables") +
  theme_bw() +
  theme(axis.text.y = element_text(color="black", size = 10),
        axis.text.x = element_text(color="black"),
        axis.title.y = element_blank(),
        axis.title.x = element_text(size = 9),
        legend.position = "none",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()) +
  facet_grid(. ~ Pais)


# Plotting results for Socio-demographic variables
sociodem_III <- regIII_df %>%
  mutate(term = factor(term, level=c("Color4", "age", "IncomeDo not know/ Prefer to not answer", "Income10 minimum wages or more", 
                                     "Income5 - 10 minimum wages","Income3 - 5 minimum wages", "Income2 - 3 minimum wages", "Income1 - 2 minimum wages", 
                                     "ReligionOthers/No Relig.", "ReligionEvangelical Traditional", "ReligionEvangelical Pentecostal or other evangelical", 
                                     "ReligionCatholic","Escolaridade_3Undergraduate or more", "Escolaridade_3High school or equivalent", "FemaleYes", "pi_cs", "pi_lf"))) %>%
  filter(term %in% demog) %>%
  mutate(term = dplyr::recode(term,
                       "pi_lf" = "Political ideology: Left", 
                       "pi_cs" = "Political ideology: Progressive", 
                       "Color4" = "Race: Black", 
                       "age" = "Age (Years)", 
                       "IncomeDo not know/ Prefer to not answer" = "    Income: Do not know / Prefer to not answer", 
                       "Income10 minimum wages or more" = "Income: 10 minimum wages or more", 
                       "Income5 - 10 minimum wages" = "Income: 5 - 10 minimum wages",
                       "Income3 - 5 minimum wages" = "Income: 3 - 5 minimum wages",
                       "Income2 - 3 minimum wages" = "Income: 2 - 3 minimum wages", 
                       "Income1 - 2 minimum wages" = "Income: 1 - 2 minimum wages", 
                       "ReligionOthers/No Relig." = "Religion: Others/No Relig.", 
                       "ReligionEvangelical Traditional" = "Religion: Evangelical Traditional", 
                       "ReligionEvangelical Pentecostal or other evangelical" = "Religion: Evangelical Pentecostal or other", 
                       "ReligionCatholic" = "Religion: Catholic",
                       "Escolaridade_3Undergraduate or more" = "Education: Undergraduate or more", 
                       "Escolaridade_3High school or equivalent" = "Education: High school or equivalent", 
                       "FemaleYes" = "Female")) %>%
  ggplot(aes(x = estimate, y = term)) +
  geom_vline(xintercept = 0, colour = gray(1/2), lty = 2, size = 0.25) +
  geom_linerange(aes(x = estimate, 
                     xmin = conf.low,
                     xmax = conf.high),
                 lwd = 0.65) +
  geom_point(aes(x = estimate, 
                 y = term),
             shape = 21,
             size = 2,
             fill = "white") +
  xlab("Estimate") +
  ggtitle("Political ideology and \n Socio-Demographic variables") +
  theme_bw() +
  theme(axis.text.y = element_text(color="black", size = 10),
        axis.text.x = element_text(color="black"),
        axis.title.y = element_blank(),
        axis.title.x = element_text(size = 9),
        legend.position = "none",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()) +
  facet_grid(. ~ Pais)

# Joining both plots
grid.arrange(indepen_III, sociodem_III,
             ncol = 1)

#------------------------------------------------------------------------------#
### Plot (IV) : Belief in Climate Change Index per Country (OLS) #####################
#------------------------------------------------------------------------------#

regIV <- feols(index_std ~ pi_lf + pi_cs + sk + ok + sc + ts + nep + ii + ei + pe +
                 Female + Escolaridade_3 + Religion + Income + age + Color4, 
               data=Base, weights= ~ ponde2, split = ~ Pais, vcov = "hetero")

names(regIV) <- countries

# Saving data frame with estimates and ci of model for each country

# Creating function to add CI to tidy function
tidy_ci <- function(x){
  tidy(x, conf.int = TRUE, conf.level = 0.95)
} 

regIV_df <- map_dfr(regIV, tidy_ci, .id = "Pais") %>%
  filter(term != "(Intercept)")

# Plotting results for Independent variables
indepen_IV <- regIV_df %>%
  mutate(term = factor(term, level=c("pe", "ei", "ii", "nep", "ts", "sc", "ok", "sk"))) %>%
  filter(term %in% indepen) %>%
  mutate(term = dplyr::recode(term,
                       "sk" = "Subjective knowledge", 
                       "ok" = "Objective knowledge", 
                       "sc" = "Scientific consensus",
                       "ts" = "Trust in scientists", 
                       "nep" = "The New Ecological Paradigm (NEP)", 
                       "ii" = "Individualism worldview", 
                       "ei" = "Egalitarianism worldview", 
                       "pe" = "Personal experience (extreme weather events)")) %>%
  ggplot(aes(x = estimate, y = term)) +
  geom_vline(xintercept = 0, colour = gray(1/2), lty = 2, linewidth = 0.25) +
  geom_linerange(aes(x = estimate, 
                     xmin = conf.low,
                     xmax = conf.high),
                 lwd = 0.65) +
  geom_point(aes(x = estimate, 
                 y = term),
             shape = 21,
             size = 2,
             fill = "white") +
  xlab("Estimate") +
  ggtitle("Psychological variables") +
  theme_bw() +
  theme(axis.text.y = element_text(color="black", size = 10),
        axis.text.x = element_text(color="black"),
        axis.title.y = element_blank(),
        axis.title.x = element_text(size = 9),
        legend.position = "none",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()) +
  facet_grid(. ~ Pais)


# Plotting results for Socio-demographic variables
sociodem_IV <- regIV_df %>%
  mutate(term = factor(term, level=c("Color4", "age", "IncomeDo not know/ Prefer to not answer", "Income10 minimum wages or more", 
                                     "Income5 - 10 minimum wages","Income3 - 5 minimum wages", "Income2 - 3 minimum wages", "Income1 - 2 minimum wages", 
                                     "ReligionOthers/No Relig.", "ReligionEvangelical Traditional", "ReligionEvangelical Pentecostal or other evangelical", 
                                     "ReligionCatholic", "Escolaridade_3Undergraduate or more", "Escolaridade_3High school or equivalent", "FemaleYes", "pi_cs", "pi_lf"))) %>%
  filter(term %in% demog) %>%
  mutate(term = dplyr::recode(term,
                       "pi_lf" = "Political ideology: Left", 
                       "pi_cs" = "Political ideology: Progressive", 
                       "Color4" = "Race: Black", 
                       "age" = "Age (Years)", 
                       "IncomeDo not know/ Prefer to not answer" = "    Income: Do not know / Prefer to not answer", 
                       "Income10 minimum wages or more" = "Income: 10 minimum wages or more", 
                       "Income5 - 10 minimum wages" = "Income: 5 - 10 minimum wages",
                       "Income3 - 5 minimum wages" = "Income: 3 - 5 minimum wages",
                       "Income2 - 3 minimum wages" = "Income: 2 - 3 minimum wages", 
                       "Income1 - 2 minimum wages" = "Income: 1 - 2 minimum wages", 
                       "ReligionOthers/No Relig." = "Religion: Others/No Relig.", 
                       "ReligionEvangelical Traditional" = "Religion: Evangelical Traditional", 
                       "ReligionEvangelical Pentecostal or other evangelical" = "Religion: Evangelical Pentecostal or other", 
                       "ReligionCatholic" = "Religion: Catholic",
                       "Escolaridade_3Undergraduate or more" = "Education: Undergraduate or more", 
                       "Escolaridade_3High school or equivalent" = "Education: High school or equivalent", 
                       "FemaleYes" = "Female")) %>%
  ggplot(aes(x = estimate, y = term)) +
  geom_vline(xintercept = 0, colour = gray(1/2), lty = 2, size = 0.25) +
  geom_linerange(aes(x = estimate, 
                     xmin = conf.low,
                     xmax = conf.high),
                 lwd = 0.65) +
  geom_point(aes(x = estimate, 
                 y = term),
             shape = 21,
             size = 2,
             fill = "white") +
  xlab("Estimate") +
  ggtitle("Political ideology and \n Socio-Demographic variables") +
  theme_bw() +
  theme(axis.text.y = element_text(color="black", size = 10),
        axis.text.x = element_text(color="black"),
        axis.title.y = element_blank(),
        axis.title.x = element_text(size = 9),
        legend.position = "none",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()) +
  facet_grid(. ~ Pais)

# Joining both plots
grid.arrange(indepen_IV, sociodem_IV,
             ncol = 1)

